library(ggplot2)
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/bayes_BF/")
sel1=list.files(pattern = "bayes_factor.txt")

group1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M28_M27","M30_M29","M26_M25","M35_M36","M18_M17","M20_M19","M22_M21","M40_M39","M50_M49")
group2=c(rep(c("DC"),5),rep("CC",4),rep("NC",4),rep(c("DC"),1))
sample_info=data.frame(group1,group2)

rt=data.frame(matrix(NA,1,4,byrow = T))
names(rt)=c("unitID","BayesFactor","group1","group2")
rt=rt[-1,]

for(i in 1:length(sel1)){
file=read.table(sel1[i],head=T,sep="\t")
file$unitID=paste(file$chrom,file$position,sep=":")
file$group1=gsub(".bayes_factor.txt","",sel1[i])
file=merge(file,sample_info,by = "group1")

tmp=data.frame(unitID=file$unitID,BayesFactor=file$BayesFactor,group1=file$group1,group2=file$group2)
rt=rbind(rt,tmp)
}
OutVals = boxplot(rt$BayesFactor)$out
rt1=rt[!rt$BayesFactor %in% OutVals,]
rt2=rt[rt$BayesFactor<50&rt$BayesFactor>10,]

factor1=factor(rt2$group1,levels = c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49","M28_M27","M30_M29","M26_M25","M35_M36","M18_M17","M20_M19","M22_M21","M40_M39"))
factor2=factor(rt2$group2,levels = c("DC","CC","NC"))

ggplot(rt2, aes(x=factor2, y=BayesFactor),color=group2) +
  geom_boxplot() +
  theme(axis.text.x=element_text(angle=50,hjust=0.5, vjust=0.5)) +
  theme(legend.position="none")

ggplot(rt2, aes(x=factor2, y=BayesFactor),color=group2) +
  geom_violin(aes(fill=factor(group2))) +
  theme(axis.text.x=element_text(angle=50,hjust=0.5, vjust=0.5)) +
  theme(legend.position="none")

test=rt[rt$group2=="NC",]
table(test$BayesFactor>10)
 